In [6]:
%matplotlib inline
Ce notebook est la traduction française du cours sur SymPy disponible entre autre sur Wakari avec quelques modifications et compléments notamment pour la résolution d'équations différentielles. Il a pour but de permettre aux étudiants de différents niveaux d'expérimenter des notions mathématiques en leur fournissant une base de code qu'ils peuvent modifier.
SymPy - est un module Python qui peut être utilisé dans un programme Python ou dans une session IPython. Il fournit de puissantes fonctionnalités de calcul symbolique.
Pour commencer à utiliser SymPy dans un programme ou un notebook Python, importer le module sympy:
In [4]:
from sympy import *
Pour obtenir des sorties mathématiques formatées $\LaTeX$ :
In [8]:
from sympy import init_printing
init_printing(use_latex=True)
In [9]:
x = Symbol('x')
In [10]:
(pi + x)**2
Out[10]:
In [11]:
# manière alternative de définir plusieurs symboles en une seule instruction
a, b, c = symbols("a, b, c")
On peut ajouter des contraintes sur les symboles lors de leur création :
In [12]:
x = Symbol('x', real=True)
In [13]:
x.is_imaginary
Out[13]:
In [6]:
x = Symbol('x', positive=True)
In [15]:
x > 0
Out[15]:
In [16]:
1+1*I
Out[16]:
In [17]:
I**2
Out[17]:
In [18]:
(1 + x * I)**2
Out[18]:
In [19]:
r1 = Rational(4,5)
r2 = Rational(5,4)
In [20]:
r1
Out[20]:
In [21]:
r1+r2
Out[21]:
In [22]:
r1/r2
Out[22]:
SymPy permet une précision arbitraire des évaluations numériques et fournit des expressions pour quelques constantes comme : pi, E, oo pour l'infini.
Pour évaluer numériquement une expression nous utilisons la fonction evalf (ou N). Elle prend un argument n qui spécifie le nombre de chiffres significatifs.
In [23]:
pi.evalf(n=50)
Out[23]:
In [24]:
E.evalf(n=4)
Out[24]:
In [25]:
y = (x + pi)**2
In [26]:
N(y, 5) # raccourci pour evalf
Out[26]:
Quand on évalue des expressions algébriques on souhaite souvent substituer un symbole par une valeur numérique. Dans SymPy cela s'effectue par la fonction subs :
In [27]:
y.subs(x, 1.5)
Out[27]:
In [28]:
N(y.subs(x, 1.5))
Out[28]:
La fonction subs permet de substituer aussi des symboles et des expressions :
In [29]:
y.subs(x, a+pi)
Out[29]:
On peut aussi combiner l'évolution d'expressions avec les tableaux de NumPy (pour tracer une fonction par ex) :
In [2]:
import numpy
In [7]:
x_vec = numpy.arange(0, 10, 0.1)
y_vec = numpy.array([N(((x + pi)**2).subs(x, xx)) for xx in x_vec])
Out[7]:
In [32]:
import matplotlib.pyplot as plt
In [33]:
fig, ax = plt.subplots()
ax.plot(x_vec, y_vec);
Une des principales utilisations d'un système de calcul symbolique est d'effectuer des manipulations algébriques d'expression. Il est possible de développer un produit ou bien de factoriser une expression. Les fonctions pour réaliser ces opérations de bases figurent dans les exemples des sections suivantes.
In [34]:
(x+1)*(x+2)*(x+3)
Out[34]:
In [35]:
expand((x+1)*(x+2)*(x+3))
Out[35]:
La fonction expand (développer) prend des mots clés en arguments pour indiquer le type de développement à réaliser. Par exemple pour développer une expression trigonomètrique on utilise l'argument trig=True :
In [36]:
sin(a+b)
Out[36]:
In [37]:
expand(sin(a+b), trig=True)
Out[37]:
In [38]:
sin(a+b)**3
Out[38]:
In [39]:
expand(sin(a+b)**3, trig=True)
Out[39]:
Lancer help(expand) pour une explication détaillée des différents types de développements disponibles.
L'opération opposée au développement est bien sur la factorisation qui s'effectue grâce à la fonction factor :
In [40]:
factor(x**3 + 6 * x**2 + 11*x + 6)
Out[40]:
In [41]:
x1, x2 = symbols("x1, x2")
In [42]:
factor(x1**2*x2 + 3*x1*x2 + x1*x2**2)
Out[42]:
In [43]:
# simplify expands a product
simplify((x+1)*(x+2)*(x+3))
Out[43]:
In [44]:
# simplify uses trigonometric identities
simplify(sin(a)**2 + cos(a)**2)
Out[44]:
In [45]:
simplify(cos(x)/sin(x))
Out[45]:
simplify permet aussi de tester l'égalité d'expressions :
In [46]:
exp1 = sin(a+b)**3
exp2 = sin(a)**3*cos(b)**3 + 3*sin(a)**2*sin(b)*cos(a)*cos(b)**2 + 3*sin(a)*sin(b)**2*cos(a)**2*cos(b) + sin(b)**3*cos(a)**3
simplify(exp1 - exp2)
Out[46]:
In [47]:
if simplify(exp1 - exp2) == 0:
print "{0} = {1}".format(exp1, exp2)
else:
print "exp1 et exp2 sont différentes"
In [48]:
f1 = 1/((a+1)*(a+2))
In [49]:
f1
Out[49]:
In [50]:
apart(f1)
Out[50]:
In [51]:
f2 = 1/(a+2) + 1/(a+3)
In [52]:
f2
Out[52]:
In [53]:
together(f2)
Out[53]:
Simplify combine les fractions mais ne factorise pas :
In [54]:
simplify(f2)
Out[54]:
In [55]:
y
Out[55]:
Dérivée première
In [56]:
diff(y**2, x)
Out[56]:
Pour des dérivées d'ordre supérieur :
In [57]:
diff(y**2, x, x) # dérivée seconde
Out[57]:
In [58]:
diff(y**2, x, 2) # dérivée seconde avec une autre syntaxe
Out[58]:
Pour calculer la dérivée d'une expression à plusieurs variables :
In [59]:
x, y, z = symbols("x,y,z")
In [60]:
f = sin(x*y) + cos(y*z)
$\frac{d^3f}{dxdy^2}$
In [61]:
diff(f, x, 1, y, 2)
Out[61]:
In [62]:
f
Out[62]:
In [63]:
integrate(f, x)
Out[63]:
En fournissant des limites pour la variable d'intégration on peut évaluer des intégrales définies :
In [64]:
integrate(f, (x, -1, 1))
Out[64]:
et aussi des intégrales impropres pour lesquelles on ne connait pas de primitive
In [65]:
x_i = numpy.arange(-5, 5, 0.1)
y_i = numpy.array([N((exp(-x**2)).subs(x, xx)) for xx in x_i])
fig2, ax2 = plt.subplots()
ax2.plot(x_i, y_i)
ax2.set_title("$e^{-x^2}$")
Out[65]:
In [66]:
integrate(exp(-x**2), (x, -oo, oo))
Out[66]:
Rappel, oo est la notation SymPy pour l'infini.
In [67]:
n = Symbol("n")
In [68]:
Sum(1/n**2, (n, 1, 10))
Out[68]:
In [69]:
Sum(1/n**2, (n,1, 10)).evalf()
Out[69]:
In [70]:
Sum(1/n**2, (n, 1, oo)).evalf()
Out[70]:
In [71]:
N(pi**2/6) # fonction zeta(2) de Riemann
Out[71]:
Les produits sont calculés de manière très semblables :
In [72]:
Product(n, (n, 1, 10)) # 10!
Out[72]:
In [73]:
limit(sin(x)/x, x, 0)
Out[73]:
On peut changer la direction d'approche du point limite par l'argument du mot clé dir :
In [74]:
limit(1/x, x, 0, dir="+")
Out[74]:
In [75]:
limit(1/x, x, 0, dir="-")
Out[75]:
In [76]:
series(exp(x), x)
Out[76]:
Par défaut le développement de l'expression s'effectue au voisinage de $x=0$, mais on peut développer la série au voisinage de toute autre valeur de $x$ en incluant explicitement cette valeur lors de l'appel à la fonction :
In [77]:
series(exp(x), x, 1)
Out[77]:
Et on peut explicitement définir jusqu'à quel ordre le développement doit être mené :
In [78]:
series(exp(x), x, 1, 10)
Out[78]:
Le développement en séries inclue l'ordre d'approximation. Ceci permet de gérer l'ordre du résultat de calculs utilisant des développements en séries d'ordres différents :
In [79]:
s1 = cos(x).series(x, 0, 5)
s1
Out[79]:
In [80]:
s2 = sin(x).series(x, 0, 2)
s2
Out[80]:
In [81]:
expand(s1 * s2)
Out[81]:
Si on ne souhaite pas afficher l'ordre on utilise la méthode removeO :
In [82]:
expand(s1.removeO() * s2.removeO())
Out[82]:
Mais cela conduit à des résultats incorrects pour des calculs avec plusieurs développements :
In [83]:
(cos(x)*sin(x)).series(x, 0, 6)
Out[83]:
In [84]:
m11, m12, m21, m22 = symbols("m11, m12, m21, m22")
b1, b2 = symbols("b1, b2")
In [85]:
A = Matrix([[m11, m12],[m21, m22]])
A
Out[85]:
In [86]:
b = Matrix([[b1], [b2]])
b
Out[86]:
Avec les instances de la classe Matrix on peut faire les opérations algébriques classiques :
In [87]:
A**2
Out[87]:
In [88]:
A * b
Out[88]:
Et calculer les déterminants et inverses :
In [89]:
A.det()
Out[89]:
In [90]:
A.inv()
Out[90]:
In [91]:
solve(x**2 - 1, x)
Out[91]:
In [92]:
solve(x**4 - x**2 - 1, x)
Out[92]:
In [93]:
expand((x-1)*(x-2)*(x-3)*(x-4)*(x-5))
Out[93]:
In [94]:
solve(x**5 - 15*x**4 + 85*x**3 - 225*x**2 + 274*x - 120, x)
Out[94]:
Système d'équations :
In [95]:
solve([x + y - 1, x - y - 1], [x,y])
Out[95]:
En termes d'autres expressions symboliques :
In [96]:
solve([x + y - a, x - y - c], [x,y])
Out[96]:
In [97]:
from sympy import Function, dsolve, Eq, Derivative, sin, cos, symbols
from sympy.abc import x
Exemple d'équation différentielle du 2e ordre
In [98]:
f = Function('f')
dsolve(Derivative(f(x), x, x) + 9*f(x), f(x))
Out[98]:
In [99]:
dsolve(diff(f(x), x, 2) + 9*f(x), f(x), hint='default', ics={f(0):0, f(1):10})
Out[99]:
In [141]:
# Essai de récupération de la valeur de la constante C1 quand une condition initiale est fournie
eqg = Symbol("eqg")
g = Function('g')
eqg = dsolve(Derivative(g(x), x) + g(x), g(x), ics={g(2): 50})
eqg
Out[141]:
In [144]:
print "g(x) est de la forme {}".format(eqg.rhs)
In [129]:
# recherche manuelle de la valeur de c1 qui vérifie la condition initiale
c1 = Symbol("c1")
c1 = solve(Eq(c1*E**(-2),50), c1)
print c1
SymPy ne sait pas résoudre cette equation différentielle non linéaire avec $h(x)^2$ :
In [112]:
h = Function('h')
try:
dsolve(Derivative(h(x), x) + 0.001*h(x)**2 - 10, h(x))
except:
print "une erreur s'est produite"
In [110]:
from scipy.integrate import odeint
def dv_dt(vec, t, k, m, g):
z, v = vec[0], vec[1]
dz = -v
dv = -k/m*v**2 + g
return [dz, dv]
vec0 = [0, 0] # conditions initiales [altitude, vitesse]
t_si = numpy.linspace (0, 30 ,150) # de 0 à 30 s, 150 points
k = 0.1 # coefficient aérodynamique
m = 80 # masse (kg)
g = 9.81 # accélération pesanteur (m/s/s)
v_si = odeint(dv_dt, vec0, t_si, args=(k, m, g))
print "vitesse finale : {0:.1f} m/s soit {1:.0f} km/h".format(v_si[-1, 1], v_si[-1, 1] * 3.6)
fig_si, ax_si = plt.subplots()
ax_si.set_title("Vitesse en chute libre")
ax_si.set_xlabel("s")
ax_si.set_ylabel("m/s")
ax_si.plot(t_si, v_si[:,1], 'b')
Out[110]: